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SUMMARY 


This report presents an overview of the probabilistic finite 
element method (PFEM) developed by the authors and their 
colleagues in recent years. The primary focus is placed on the 
development of PFEM for both structural mechanics problems and 
fracture mechanics problems. The perturbation techniques are 
used as major tools for the analytical derivation. The report 
specifically covers the following: (1) representation and 

discretization of random fields; (2) development of PFEM for the 
general linear transient problem and nonlinear elasticity using 
Hu-Washizu variational principle; (3) computational aspects; (4) 
discussions of the application of PFEM to the reliability 
analysis of both brittle fracture and fatigue; and (5) a novel 
stochastic computational tool based on stochastic boundary 
element (SBEM) . Results are obtained for the reliability index 
and corresponding probability of failure for: (1) fatigue crack 

growth; (2) defect geometry; (3) fatigue parameters and; (4) 
applied loads. These results show that initial defect is a 
critical parameter. 
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1. INTRODUCTION 

It is becoming increasingly evident that traditional and especially deterministic methods will 
„oc be sufLientTproSy design advanced structures or structural components subjected to a 

of complex Ind cyclic toading conditions. Due to uncertainty in loading condition , 
material behavior geometric configuration, and supports, the stochastic computational mechanics, 
XSccount“feTtf.ese uncertain aspects, has to be applied to provide rauonal teltabth y 
anaNs^ C< tmdm describe the behavior of the structure. The fundamentals of stochastic 
computational mechanics and its application to the analysis of uncertain structural systems are 
summarized and recapitulated in the book (Liu and Belytschko, 1989). 

While the theory of statistics and structural reliability has been used successfully in 
modeling the uncertain nature of structures, load environments, and in computing die pro a 1 ty 
failure its application is only limited in simple structures with linear constitutive Savior ® . 

the comSy in the geometry, external loads, and nonlinear matenal behavior, more advanced 
?ompuTti3 tools such as Finite Element Methods (FEMs) or Boundaiy Integral Equation 
Methods (BIEMs), have to be employed to provide the necessary computational framew 
analyzing structural response. The combination of these advanced computational tools with the 
theory ©^statistics and structural reliability has become a rational way for the safety assessment and 
nSainw. characterization of complex structures. In this Chapter, attenuation is focused on the 
development of Probabilistic Finite Element Method (PFEM), which combines the finiteelemen 
method with statistics and reliability methods, and its application to linear, nonlinear structural 
mechanics problems and fracture mechanics problems. The novel computational tool based on the 
Stochastic Boundary Element Method is also given for the reliability analysis of a curvilinear 

fatigue crack growth. 

The existing PFEMs have been applied to solve for two types of problems: (1) 
determination of the response uncertainty in terms of the means, vanance and correlation 
coefficients- (2) determination the probability of failure associated with prescribed limit states. 

rScSd SSSstic moments of a 'response are no. sufficient for a coopte ntobtoy 
analysis, these moments offer useful statistical information and sente as amcasures of reUabilit^ 
Furthermore, due to the lack of multivariate distribution function of random variables, a 
meaningful risk assessment and failure analysis may not be feasible. 

The perturbation method has been used extensively in developing PFEM due to ns 
simplicity and versatility. Cambou (1975) appears to have been the first to apply the first order 
nerturbation method for the finite element solution of linear static problems with loading and 
SmSS Baecher and Ingra (1981) also used the same techniques for settlements 
predictions The perturbation method in conjunction with finite element method has also been 
adopted by Handaand Anderson (1981) for static problems of beana and frame smictures, by Ishn 
and Suzuki (1987) for slope stability reliability analysis, and by Righetn and Harrop-Williams 
(1988) for static stress analysis for soils. The accuracy, convergence and computational efficiency 
of the perturbation method have been compared with those from Neumann expansion method and 
direct Monte Carlo Simulation (MCS) method (Shinozuka and Yamazaki, 1988, “ 

Deodatis 1988) The PFEM based on the second-order perturbation approximation has been 
introduced by Hisada and Nakagiri (1981 and 1985) for static problems and for eigenvalue 

problems. 

Extensive research on the PFEM has been developed by the authors and their coUeagues at 
Northwestern University in the recent years. The PFEM based on the second-order l^rturbation has 
been developed to estimate the statistic moments of the response for linear static problems (Liu et al., 
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1986a), nonlinear dynamic problems (Liu et al., 1986b), and inelastic problems (Liu et al., 1987). 
The .formulation based on the single-field variational principle has been extended by Liu et al. 
(1988a) to the three-field Hu-Washizu variational principle formulation, which has far greater 
versatility. The numerical instability resulting from the secular terms in the perturbation has been 
removed by Liu et al. (1988b) based on Fourrier analysis. The perturbation methods have been 
shown to provide efficient and accurate results for small random fluctuations in the ra ndo m 
parameters. An extensive review on the application of perturbation methods in developing PFEM 
has been given by Benaroya and Rehak (1988). 

The finite element method coupled with the First and Second-Order Reliability Methods 
(FORM and SORM) has been developed by Der Kiureghian and Ke (1985, 1988) for linear 
structural problems and Liu and Der Kiureghian (1991) for geometrically non-linear problems. The 
most critical step in this method is the development of an efficient search algorithm for locating the 
point at which the response surface is to be expanded in a first or second order Taylor senes. This 
point is obtained by an iterative optimization algorithm, which involves repeated computation of the 
limit state function and response derivatives. Unlike the method of direct differentiation (Der 
Kiureghian and Ke, 1988; Liu and Der Kiureghian, 1991; Zhang and Der Kiureghian, 1991), .the 
PFEM based on the perturbation approximation in conjunction with FORM has been developed by 
Besterfield et al (1990, 1991) for the reliability analysis of brittle fracture and fatigue. In a slightly 
different context, the PFEM has been developed by Faravelli (1986, 1989) that couples response 
surface approach with deterministic finite element formulation. The finite element simulation 
coupled with the polynomial response surface fitting has been also proposed by Grigoriu (1982). 
Using a deterministic finite element code and finite differences, an advanced algorithm based on the 
Fast Probability Integration (FPI) has been developed by Wu et al (1990) to generate the entire part 
of the Cumulative Distribution Function (CDF) of the response. The performance of the FPI based 
on either advanced mean-value method or advanced mean value first-order method has been 
demonstrated by Cruse et al. (1988) through the reliability analysis of turbine blade. 

In addition to the PFEM, the Stochastic Boundary Element Method (SBEM) has been 
developed and adopted recently by researchers. The SBEM that combines the deterministic 
boundary element method with perturbation expansions has been developed by Ettouney et al. 
(1989) and Dasgupta (1992) for the determination of the statistic moments of both displacements and 
tractions. Most recently, the authors have developed the SBEM, which combines the mixed 
boundary integral equation method (Lua et al., 1992c) with FORM, for the study of probabilistic 
fatigue crack growth (Lua et al., 1992d). 

This Chapter presents an overview of the PFEM developed by the authors and their 
colleagues in the recent years. The primary focus is placed on the development of PFEM for both 
structural mechanics problems and fracture mechanics problems. The perturbation techniques are 
used as major tools for the analytical derivation. The remainder of this Chapter is organized as 
follows. In Section 2, the representation and discretization of random fields are presented. The 
development of PFEM for the general linear transient problem and nonlinear elasticity using Hu- 
Washizu variational principle are given in Section 3, and 4, respectively. The computational aspects 
are discussed in Section 5. The application of PFEM to the reliability analysis of both brittle fracture 
and fatigue is given in Section 6. A novel stochastic computational tool based on SBEM is 
presented in Section 7. The final conclusions are drawn in Section 8. 

2. RANDOM FIELD DISCRETIZATION 


2.1 Background 

The randomness of a stochastic system can be described in three forms, random variables, 
random process in space, and random process in time. The random process in space is also called 
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random field. The aspects of random fields and its application to engineering problems are given 
by Vanmarcke (1984). Various methods have been used to the numerical representation ot random 
processes. The statistical characterization for thunderstorm winds has been given by Twisdale and 
Dunn (1983) and Twisdale and Vickery (1991). The spectral representation of random processes 
based on computer simulation has been proposed by Shinozuka (1987). 

The spatial variability of mechanical properties of a system and the intensity of a distributed 
load can conveniently be represented by means of random fields. Due to the discrete nature of the 
finite element formulation, the random field must also be discretized into random variables. This 
process is commonly known as random field discretization. Various methods have been developed 
in the representation of random fields. They are: the midpoint method (Hisada and Nakagin, 1985; 
Der Kiureghian and Ke, 1988; Yamazaki et al., 1988), the spatial averaging method (Vanmarcke 
and Grigoriu, 1983), and series expansion method (Lawrence, 1987; Spanos and Ghanem, 1988). 
In this section, the interpolation method (Liu et al., 1986a) is described. In this method, the 
random field is represented by a set of deterministic shape functions and the random nodal values 
of the field. The size of the random field element is controlled by the correlation length the fieid 
and the stability of the probability transformation used in the reliability methods (FORM and 
SORM) The random field mesh should be so fine to capture the fluctuation of the random field. 
On the other hand, the random field mesh should not be so small that high correlated stochastic 
variables of adjacent elements cause numerical instability in the probability transformanon, which 
is required in the reliability methods (FORM and SORM). As suggested by Der Kiureghian 
(1985), two separate meshes for the finite element and for random fields have to be used in the 
numerical implementation. 

Since the computational effect in the determination of response derivatives or sensitivities is 
proportional to the number of random variables, it is desirable to use as few random variables as 
possible to represent a random field. To achieve this goal, the transformation of the original 
random variables into a set of uncorrelated random variables has been introduced by Liu et al. 
(1986a) through an eigenvalue orthogonalization procedure. Comparison with a Monte-Carlo 
simulation demonstrates that a few of these uncorrelated variables with larger eigenvalues are 
sufficient for the accurate representation of the random field. This technique along with other 
computational aspects is presented in Section 5. 

2.2 Interpolation Method 

Let b(x) represent the random field. In PFEM, b(x) is approximated by 


b(x) = N i (x) b i 


where N.(x) represent the shape functions and b. the discretized values of b(x) at x., i — 1, ..., q. 
It follows from Eq. (2-1) that 

N.(x) db. ( 2 - 2 ) 

i i 



and 
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db 2 (x) 


• i 


i j = 1 


N.(x) N.(x) db. db. 
J 1 J 


(2-3) 


where 


db. = b. - b. (2-4) 

l l i 


and b. represent the mean values of b. (also denoted by the expectation operator E[.]). From Eq. 
(2-1) the expectation and the covariance of b(x) are, by definition. 


+oo 


E[b(x)]= J b(x) f(b) db 


(2-5) 


= i n. 

i = 1 1 


(x) E[b.] 


( 2 - 6 ) 


and 


+oo 


Cov[b(x k ) , b(x m )]= 


[b(x fc ) - b(x fc )] [b(x m ) - b(x m )] f(b) db 


(2-7) 


- i 


i j = 1 


N.(x. ) N.(x ) Cov(b. , b.) 
i k j m l j 


( 2 - 8 ) 


respectively, where f(b) is the multivariate probability density function; x^ and x m are any two 
points in the domain of x. 

From second-moment analysis, the mean of any function S[b(x),x] at any point x^ and the 
covariance of the function between any two points x^ and x^ can be written as 

^Cov^.b) (2-9) 

i j = 1 1 J 

and 


5 


( 2 - 10 ) 


Cov ( s k • S J “ 


_ A 


<L 9S. 

I w 


( - \ 

as 

m 


3b. 


i“iv iA u 


Cov(b. , b.) 


where 


S k = S[b(x),x k ] 


( 2 - 11 ) 


and the superposed bar implies evaluaoon at b. The error in Eqs. (2-9) and (2-iO) anses from. ( ) 
the mmcationof higher order moments and (ii) the discretization of the random ficldb(x) i by the 
finite vector b. If the randomness in b(x) is small, then the first error will be small for a smooth 
function and the second-moment analysis is applicable. The error due to discretization m Eqs. ( - 
6) and (2-8) has been studied by Liu et al. (1987) . 


When the random field discretization is coupled with a FEM discretization, as inPFEM, q 
need not be equal to the number of finite elements NUMEL and the shape function ^(x) need not 
be the same as the finite element interpolants for the displacement field As indicated before, two 
meshes, one depending on structural topology, the other on correlation length can be employed to 
im prove the computational efficiency. 

3. PROBABILISTIC FINITE ELEMENTS FOR LINEAR PROBLEMS 

The probabilistic finite element method (PFEM) is used to study of systems with paramemc 
uncertainties in both the unknown function and mathematical operators acting on loads can 

be either deterministic or random. In this section, the second-order perturbation is employed to 
develop PFEM for a general linear transient problem. By applying the second-order perturbation, 
the random linear system equations can be replaced by a finite number of identical detenmnistic 
system equations up to second-order. The effective load in each of these equations depends on the 
randomness of the system and the solutions of the all the lower order equations. 

Due to the space limitation, the review on the deterministic finite element method is not 
given here. The state-of-the-art of finite element techniques can be found in the review article by 
Noor (1991). Using either the single field variational principle or Galerkin formulation, the 
discretized linear equations of motion are 


M a(b, t) + K(b) d(b, t) = f(b, t) 


(3-1) 


where M and K(b) are the (neq x neq) global mass and stiffness matrices, respectively; a(b, t), 

d(b, t), and f(b, t) are the (neq x 1) nodal acceleration, displacement, and force vectors, 
respectively; neq is the number of equations; and b is a q-dimensional discretized random variable 
vector, i.e., b. = b(x.), where x. is the spatial coordinate vector. The mass is usually assumed to 

be deterministic whereas the probabilistic distributions for the stiffness and external force are 
represented by a generalized covariance matrix, Cov(b., b.), i, j = 1, .... q. It is worth noting that 

the stiffness matrix can be expressed in terms of the generalized gradient matnx, B(x), and the 
material response matrix, D(b, x). In this formulation, the random vector, b, can represent a 
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random material property (e.g., Young's modulus) and/or a random load. 

The application of second-moment analysis to develop PFEM involves expanding all 

random functions about the mean value of the random vector b, denoted by b, via Taylor series 
and retaining only terms up to second-order terms. That is, for a small parameter, the random 

displacement function d(b, t) is expanded about b via a second-order perturbation as follows: 

d(b, t) = d(t) + ^ d (t)Ab + \ £ ^ d ” b .b. (t) Ab i Ab j (3 ‘ 2) 

i=l i i j=l i j 

where d(t), d (t), and d fe b (t) represent the mean displacement, the first-order variation of 

i i j 

displacement with respect to b. evaluated at b, and the second-order variation of displacement with 

respect to b. and b. evaluated at b, respectively and Ab. represents the first-order variation of b. 

about b.. In a similar manner, K(b), a(b, t), and f(b, t) are also expanded about b. via a second- 
1 1 

order perturbation. Substitution of the second-order perturbations of the random function d(b, t), 

K(b), a(b, t) and f(b, t) into (3-1) and collecting terms of order 1, £, and % yields the following 
equations: 

Zeroth-Order Equations 

M a(t) + Kd(t) = f(t) (3-3) 

First-Order Equations (for each Ab., i = 1, .... q) 

Ma b (t) + Kd b (t) = F b (d,t) (3-4) 

l i i 

where 

F b (d, t) = f b (t) - R b d(t) (3-5) 

i i i 

Second-Order Equations (i and j are summed from 1 to q) 
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(3-6) 


M a 2 (t) + K d 2 (t) = F 2 (d, t) 


where 


F 2 (d, t) 


1 7 


^ t 2 f b.b. 

i j=l ij 


(t) 


i J 


K b d b (t) } Cov(b., b.) 
i j 


a 2 (t) 


-it 


2 t* a b.b. ft) Cov(b i’ V 

i j=l 1 J 


d 2 (t) 


i 


5 L d b.b w Cov(b i- V 

i j=l 1 J 


(3-7) 


(3-8) 


(3-9) 


The solution process for Eqs. (3-3) through (3-9) can be performed in parallel since only one 
effective stiffness matrix needs to be formed. Therefore, the total solution requires one 
factorization of the effective stiffness matrix and q+2 forward reductions and back substitutions ot 

an (neq x neq) system of linear equations to obtain the zeroth-, first-, and second-order solutions. 

To illustrate the performance of PFEM, a simple two degree of freedom spring-mass system 
is presented here. The computed results are compared with those obtained using (1) Monte Carlo 
Simulation (MCS) and (2) Hermite-Gauss Quadrature (HGQ) schemes. The problem is depicted 
in Fig. 7-1 . A sinusoidal vector forcing function is used: 


F(t) 


■[ 


0.0 

25.0 x 1 0 6 sin 2000 


.] 


(3-10) 


The random spring constants K x and K 2 are normally distributed with a coefficient of variation 

(i e < 7 /i i) equal to 0.05. The mean spring constants are 24 X 10 6 and 12 X 10 6 , respectively. The 
deterministic masses /n, and m 2 are 0.372 and 0.248, respectively. A stiffness-proportional 
damping of 3% is included. The probabilistic equationsderived earlier are solved by the implicit 

Newmark-P method (Ma, 1986). The mean amplitude d x is depicted in Fig. 7-2 for all the three 
numerical methods-PFEM, HGQ and MCS. The PFEM solution compares very well with the 
other two methods. For the variance of d\ the PFEM solution, plotted in Fig. 7-3, seems to 

overshoot the variance at large time. The ±3(7 bounds for the displacement d\ is plotted in Fig. 7- 
4. 

4. PROBABILISTIC FINITE ELEMENTS FOR NON-LINEAR PROBLEMS 


The probabilistic finite element method has been developed in the previous section using the 
single-filed variational principle. Due to the direct stiffness matrix approach used, it can be only 
applied to solve a limited number of problems with uncertainty in loading and material properties. 
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In order to consistently handle problems with randomness in the equilibrium equations, domain, 
and boundary conditions, the three-field Hu-Washizu variational principle will be employed to 
develop PFEM. An additional advantage of using the Hu-Washizu variational principle involves 
the elimination of the locking phenomena (Belytschko and Bachrach, 1986) and suppression of 
hourglass modes (Belytschko et al., 1984). Solution of three stationary conditions for the 
compatibility relation, constitutive law, and equilibrium yield the variations in displacement, strain 
and stress. The statistics such as expectation, autocovariance, and correlation of displacement, 
strain and stress are then determined. 

Using matrix notation, the Hu-Washizu Variational Principle (HWVP) for nonlinear 
problems adopted in this Section is (see Liu et al., 1988c) 


| §£ T (\|/ - a) dQ - J 8a (e-Vu)dQ 
Q Q 


f T 

f T f 

J 8(Vu) a dQ - 

J 8u F dQ - J 


T 

5u h dT 


= 0 


a 


Q 


3Q, 


(4-1) 


where e, a, and u are independent random field variables representing the nonsymmetric measure 
of the strain, first Piola-Kirchhoff stress, and displacement, respectively; \|i is a nonlinear function 
of the deformation gradient; and a superscript T represents the transpose. In Eq. (4-1): Q, 3Q^, 

F, h, and Vu represent the domain, traction boundary, body force vector, prescribed traction 

vector, and the nonsymmetric part of the displacement gradient, respectively; 8 represents a virtual 
quantity. The surface and volume integrals in Eq. (4-1) can be expressed via parametric 
representation: 


dr = J dA and dQ = J dR (4-2) 

s v 

respectively, where J s and J y represent the surface and volume Jacobians, respectively; and R and 

A represent the reference domain and boundary, respectively. Random domains and boundaries 
are incorporated into the formulation through randomness in the gradient operator and Jacobians. 
The application of second-order perturbation techniques in the HWVP involves the expansion of all 

random functions about the mean value of the random field b(x) denoted by b(x) and retaining 
only up to second-order terms, i.e., for a given small parameter = the scale of randomness in 

b(x)J, the random function y is expanded about b at a given point x in the reference domain as 
follows: 

V - V° + 5 (V + C°e ) + £ (V + C e + C°e") (4-3) 

where the superscripts nought, prime, and double prime represent the random functions evaluated 

at b, the first-order variation due to variations in b, and the second-order variation, respectively. 
The first elasticity tensor, C in Eq. (4-3) is given by 


9 



(4-4) 


d 2 ™ 

% m - sg^; 

where W is the strain energy density function; and G is the deformation gradient. Similarly, the 

T 

rest of random functions e, a, F, h, J s , J v , Vu and 8(Vu) can also be expressed as second 
order perturbations (see Liu et al„ 1988a). After substituting the second order perturbations of all 
these random functions into Eq.(4-1), the following three equations for the zeroth-, first-, and 
second-order nonlinear PHWVP are obtained: 

Zeroth-Order Variational Principle 


J 8e T (/ - a 0 ) dR - J 6a (e° - V°u°) dR 

R R 

+ J 5(V°u) T a° dR - J 5u F° J y dR - J 5u h w J"dA =0 


R R 

First-Order Variational Principle (£ terms) 




T . 0 ,o 

s 


(4-5) 


1 


5s T [(C°e - a ) J® + v J® + (¥° - o°) J y ] dR 


J 


5a T [(e - V°u - v’u°) J° + (e° - v°u°) J ] dR 


C n T > ’Ton o T 0 

+ J {(6(V°u) a + 6(V u) a U ]r + 5(V u u) a u J v )dR 


- f 8 u T (F J° + F° J ) dR - f6u (h + h° J )dA = 0 
r v V 9R, 


(4-6) 


Second-Order Variational Principle (^ terms) 


/ 


X A «« it A t t M 0 O' tit 0 0 11 

8e [(Ce -a )r + (C£+\|t)J v + (C £+\|/-o)J v + (y -a)J v ]dR 
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- J 5<j T [(£ - V°u - V u - V u°) J y + (e - V°u - V u°) J y + (e° - V°u°) J y ] dR 
R 

C n T - > T -Too 

+ { [8(V u) a + 5(V u) a + 5(V u) a ] J y 

R 

n T * ' T o 1 o T n " 

+ [8(Vti) a + 8(V u) a] J +8(Vu) a J } dR 


•J 


T " a 11 0 " 

5u (F J + F J + F u J ) dR 

V V V 


- fsu 7 (h" J° + h’ / + h° j")dA = 0 (4-7) 

dR, 

n 

It should be noted that all random functions with the superscript (" ) or (’ ) in Eqs. (4-6)-(4-7) are, 
in general, described through spatial expectation and autocovanance functions. Therefore, in 
addition to the usual finite element approximation of the displacement field, the random fields are 
also discretized with q shape functions. To be consistent with the finite element approximation and 

to mai n tai n the accuracy of the discretized random field [i.e., b(x)], the random functions y, C, F, 

h, J , and J , which are, in general, functions of b(x) and x, are first discretized with the same q 
s v 

shape functions as the random fields. For example, the finite element approximation of C is given 
by 


c = c°+ + ^ 2 c 


(4-8) 


or 


i 


c - 2 , ♦,<*)«:?+ $cj+ 4 2 C.') 

1=1 


(4-9) 


,0 


where; <)>j(x) are the q shape functions; Cj denotes the Ith nodal value of C evaluated at b; Cj 

denotes the first-order variation of C(Xj, b) due to variations Ab.; and Cj denotes the second-order 
variation. The last two are then expanded in terms of the random variables bj and given by 


C,= t (Cp,Abi 

i=l 


(4-10) 


11 


and 



(C ").. Ab. Ab. 
v I ij i J 


(4-11) 


respectively. The factor 1/2 is included in order to be consistent with the second-order Taylor 

series expansion. The nodal values (C^ and (C^ can be obtained by partial differentiation of C 

or by a least-square fit to the actual data. Similar definitions can be developed for the rest of 
random functions (see Liu et al., 1988a). 

Substituting the given approximation of all random functions into the zero-order, first- 
order, and second-order PHWVPs (Eqs. (4-5)-(4-7)), and using the three stationary conditions 
(strain-displacement, stress-strain, and equilibrium), the zeroth, first and second order equations 
can be obtained (see Liu et al., 1988a). The zeroth-order equations require an iterative solution 
technique but the first-order and second-order equation are linear. After determining the zeroth-, 
first-, and second- order solutions, the expectations and autocovanance matnces tor the 

displacements, strains, and stresses can be obtained. . 

The applicability and effectiveness of the PFEM for nonlinear problems was demonstrated 
by Liu et al. (1988a) through the problem of a cantilever beam subjected to large deflection. The 
Saint Venant-Kirchhoff model for nonlinear elasticity with randomness in the external force, beam 
height, and material properties were considered. The probabilistic distribution for displacement, 
strain and stress were also presented. 

To reduce computational effort, the random variables can be transformed to the uncorrelated 
normal form by an eigenvalue problem as shown below. 

5. COMPUTATIONAL ASPECTS 


5.1. Random Variable Transformation . 

The mean and covariance can be obtained from the equations in Section 4. However, the 
number of derivatives to be evaluated is proportional to q(q+l)/2, where q is the number of 
random variables. To reduce computations, an eigenvalue orthogonalization procedure, which is 
similar to the modal analysis in structural dynamics, can be employed. The full covariance matrix 
Cov(b. , b.) is transformed to a diagonal variance matrix Var (c. , c.) such that 


and 


Var (c. , c.) = 0 
i J 


Var (c. , Cj) = Var (c.) 


for i * j 


for i = j 


(5-1) 


(5-2) 


Therefore, the number of evaluations is proportional to q. The above is achieved through the 
eigenproblem: 


= T* A 


(5-3) 
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where the Q and A matrices denote Cov(b. , b.) and Var (c. , c^), respectively; is a constant q x 
q fundamental matrix with the following properties: 


T 

yvj/ 


T 

¥ V = I 


(5-4) 


T 

a = y ay 

and 


5 = ^0 or 


T 

c = V b 


(5-5) 

(5-6) 


I is the q x q identity matrix and c is the transformed q x 1 vector of random variables. Thus, the 
discretized random vector b is transformed to an uncorrelated random vector c with the variance of 

c as the eigenvalues of Q in Eq. (5-3). 

With Eqs. (5-5) and (5-6), the mixed derivatives appearing in Section 5 reduce to second 
derivatives and Var (b. , b.) reduces to Var (c.). Thus, the mean of any function S[b(x), x] at any 

point xfc and the coveriance of the function between any two points xfc and x m can be written as 


q a 2 s, 


EISJ = s k + 5 X — T VaI{C i ) 


and 


Cov[S k , SJ 


i = 1 3c. 


q fas ^ 

x y. 

i j = 1 V 1 J 


f - \ 
as 

m 


V^i J 


Var(c.) 


(5-7) 


(5-8) 


respectively. . ,. 

It is observed that for one-dimensional random fields, as the correlation length increases 
from zero to a large value, the number of largest eigenvalues n, n < q necessary to evaluate the 
mean and covariance in Eqs. (5-7) and (5-8) to a specified accuracy, decreases from q to 1. When 
the correlation length is zero the random field is uncorrelated and all q eigenvalues are dominant. 
As the field is uncorrelated, all q random variables are necessary to represent the randomness of the 
field. As the correlation length increases the number of dominant eigenvalues decreases. 
Eventually, for a very large correlation length the random field is closely correlated and there is just 
one dominant eigenvalue. As the field is closely correlated, only one random variable, 
corresponding to the largest eigenvalue, is sufficient to represent the randomness of the field. This 
feature, when present, can easily be exploited to reduce the computations. The value of n can be 
chosen based on the distribution of the eigenvalues before solving the PFEM equations. The 
eigenvalues here can be interpreted as weighting factors for the corresponding mode shapes 
necessary to represent the covariance structure; a large eigenvalue means a dominant mode and vice 
versa. Results of the eigenvalue distribution and selection of n, for beam problem and a bar 
problem, are discussed in Liu et al. (1986a, 1987). 


5.2. Adjoint Method in PFEM 

Consider a typical function Il(c, d) involving the displacements d and the random variables 
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c. Chain differentiation yields 

[nl = n + Tl. d ( 5 ' 9 ) 

1 J c. c. d c. 

i i i 

where the subscript denotes the derivative with respect to c., and 


Tl. = ( n H - n H • • ’ n H 
d d l d 2 d NEQ 


) 


Using the first-order equation of PFEM in the transformed space, i.e.. 


(5-10) 


i i 
E q. (5-9) becomes 

[n] =n + elk'V ( 5 - 12 > 

1 J c. c. d c. 
i i i 

Usually, in the direct method, the above equation is evaluatedfor each random variable c., 
involving 'n' solutions of the linear equation (5-12). In the adjoint method, X, is selected to satisfy 


k x = n < 5 - 13 ) 

d 

Then, Eq. (5-12) can be rewritten as 

[nl = n + x. T f (5- 14 ) 

L J c. c. c. 

1 1 1 

The adjoint problem, Eq. (5-13), is solved only once in this method. In the direct method, 'n' 
solutions of Eq. (5-11) are required. This is the advantage of the adjoint method over the direct 

method. Both methods require 'n' inner products with f in Eqs. (5-9) and (5-14), respectively. 

V • 

1 

However, it has been shown that when the number of functions is more than the number of 
random variables, the computational advantage of the adjoint method is lost (Liu et al., 1988d). 
By solving 'q' adjoint problems, the second order sensitivites can also be evaluated. It should be 
noted that the adjoint method is applicable to nonlinear problems as well, as the first and second 
order equations are still linear. 

5.3 Parallel Computing in PFEM 

Recent advances in computing hardware and software, have made multiprocessing in 
general and parallel processing in particular a viable and attractive technology. Parallel processing 
provides an opportunity to improve computing efficiency by orders of magnitude. Probabilistic 
computational mechanics exhibits several inherent levels of both coarse-and finite-grained 
parallelism. It is imperative to develop the computational strategies and algorithms to maximize 
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parallel processing efficiency and minimize parallel overhead. The parallelism in both the 
probabilistic computations and the structural mechanics computations has been explored by Sues 
et. al (1991,1992a). The implementation of parallel processing to general probabilistic structural 
analysis problems has been studied by Sues et al. (1991, 1992a). The parallel computing for 
probabilistic fatigue analysis has been implemented by Sues et al. (1992b) on both shared and 
distributed memory machines. 

The parallel implementation of PFEM can be easily achieved in the solution of the first- 
order equations (sensitivity analysis). As shown from Eqs. (3-3)-(3-9), only one effective 
stiffness matrix needs to be formulated. Once the zeroth-order solution is obtained, q equations 
(Eq. (3-4)) can be solved in parallel to determine the response derivatives. Multiple levels of 
parallelism can be achieved if the substructuring (Komzsik and Rose, 1991), domain 
decomposition (Chan, et al., 1989) and operator splitting (Sues and Chen, 1992c) are also 
employed in PFEM. 


6. APPLICATION OF PFEM TO THE RELIABILITY ANALYSIS OF BRITTLE FRACTURE 
AND FATIGUE 

6.1 Introduction 

In the previous section, the probabilistic finite element method, which is based on the 
second-order perturbation, has been formulated to quantify the statistical moments of the response of 
a stochastic system. In this section, the PFEM coupled with the first-order reliability method is 
developed to determine the reliability of brittle fracture and fatigue crack growth. The constrained 
optimization problem is formulated to calculate the reliability index. A Lagrange multiplier technique 
along with gradient projection algorithms is used to solve the constrained optimization problem. 


Fracture and fatigue have become important factors in the structural design and safety of 
aging structures. The failure of an aging structure is usually resulted from microdefects activation, 
propagation, and formation of major cracks. Due to the ramdomness in the configuration of 
microdefects and uncertainty in the failure mechanism, the probabilistic fracture mechanics (PFM), 
which combines the fracture mechanics with the stochastic methods, provides a useful tool to 
address problems with large uncertainty. 

The reliability analysis of flawed structures will here be classified into two groups of 
problems: 

1. brittle material problems, where the material contains flaws with the random location and 
orientation. The major question is the reliability of the structure in the presence of these flaws. 

2. ductile material problems, where failure is expected to result from the growth of a critical 
flaw until it can lead to failure of the structure. 

The first category of the problem has been addressed recently by Lua et al. (1992a, 1992b) in 
quantifying the inherent statistical distribution of the fracture toughness of a multi-phase brittle 
material. The second question is of particular relevance to the safety of aging structures and 
nondestructive evaluation techniques. Because the threshold of detection is substantially greater 
than flaws sizes which may lead to failure over the course of time, inspection cycles should be set 
so that the reliability of an aging structure remains acceptable in these circumstances. Although a 
deterministic analysis can obtain an estimate of the fatigue life, the uncertainties in crack growth 
rates and the initial crack lengths detract from the usefulness of such solutions. 
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In this section, the application of the PFEM and the first-order reliability analysis for the 
brittle fracture and fatigue is presented. A brief review on the reliability methods is given first. 
The fusion of the PFEM and reliability analysis for probabilistic fracture and fatigue are then 
presented Performance of the methodology developed is demonstrated on example problems. 

6.2 Reliability Analysis 

Reliability analysis is used to determine what is the likelihood that a smicture subjected to 
uncertain loads, material properties and geometry, will satisfy a limit state entenoii. Several text 
books and monographs on the methods and application of the reliability theory have been vmtten, 
i.e. [Ditlevsen (1981), Ang and Tang (1984), Augusti et al. (1984), Madsen et al. 0986) and 
Melchers (1987)1. An overview of the first-order and seconder-order reliability method 
(FORM/SORM) as well as various Monte Carlo simulation schemes has been given by Bjerager 
(1989) As the PFEM provides a powerful computational tool to determine first- and second- 
moment of random parameters, the second-moment reliability methods can be easily combined 
with the PFEM to obtain measures of the reliability of the structural systems. 

Throughout this section the uncertainties-in load, material properties, component geometry 
and crack geometry- are represented by a q-dimensional vector of random variables denoted by b - 

[b , ..., b ] T A random variable reliability problem is described by a performance function, g(b), 

which is a continuous measure of the ability of a system to perform as designed. Three states of a 
system, namely, the limit-state surface, the failure state, and the safe state, are defined by 


g(b)=0 


g(b)<0 and g(b)>0 


( 6 - 1 ) 


respectively. The probability of failure is given by 


r ’I 


P f = I f B (b) db (6-2) 

J iOOSO 

where fg(b) is the multivariate density function of b. Two difficulties are associated with Eq. (6- 

2). First, the domain of integration ( g(b) < 0) is an implicit function of the random vector b. 
Second, standard numerical integration of this multiple integral is prohibitively complicated as the 
number of random variables becomes large. Two approaches- MCS 
approximation methods such as the first or second order reliability method (FORM or SORM)- 
have been employed extensively to calculate Eq. (6-2). In the first-order reliability method 
(FORM) the limit-state surface in the standard normal space is represented by the tangent 
hyperplane at the design point. In the second-order reliability method (SORM), the limit-state 
surface in the standard normal space is replaced by a quadratic surface, which is tangent at the 
design point. While MCS is completely general, it is very expensive and time-consuming for small 
probabilities of failure, which is the major concern in reliability engineering. FORM and SORM 
are more accurate and efficient for extreme probability of failure (e.g., 0.0001 or 0.9999 
cumulative probability), however implementation can be more complex. In the present study, the 
FORM is applied to predict the reliability of a flawed component . 


In order to make use of the properties of the standard normal space (rotationally symmetric 
and exponential decay), a transformation is introduced to map the original random variables b to a 
set of standard, uncorrelated normal variables r. Eq. (6-2) in the r-space becomes 
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P f = I (IrcJ'^expj-i-r'rjdr 

Jg( ^ ^ 0 


(6-3) 


where ( / denotes the transpose of a vector or a matrix, and g(r) (- g(b(r))) is the performance 
function in the transformed r-space. FORM approximates the calculation in Eq.(6-3) as follows: 
first the point r* on the limit-state surface (g(r) = 0), which has the minimum distance to the origin, 
is found through an iterative algorithm, then the limit-state surface at the design point r* is 
replaced with a tangent hyperplane given by 

g(r) = (r, - r* ) (6-4) 

3ri 

The resulting first order approximation to Eq. (6-3) is 


r 

, (2rc)' q/2 exp(- ^ r l r)dr = d>(-(}) 

^ to - r')<0 2 


where the reliability index (3 is defined by 

P = (r*V) 1/2 (6-6) 

and <I>( ) is the standard normal cumulative probability. The step to determine the most probable 
point (r* )on the failure surface is the most critical in the reliability analysis. It generally requires to 
form an iteration and optimization scheme to calculate the gradients of the performance function. 

In this paper, the reliability index (3 is determined by solving the following optimization 
problem in r-space, i.e., 

p = min (r T r) 1/2 , subject to g(r) = 0 (6-7) 

The optimization can be solved using any general non-linear optimization algorithm such as HL-RF 
method (Hasofer and Lind, 1974; Hohenbichler and Rackwitz, 1981; Rackwitz and Fiessler, 
1978), gradient projection method (Haug and Arora, 1979) and the modified HL-RF method (Der 
Kiureghian and Liu, 1988). A fast convergence rate is essential in selecting an iteration method. 

The second order reliability method based on the second order Taylor expansion of the 
failure surface is given by Fiessler et al. (1979), Breitung (1984), Der Kiureghian et al. (1987), 
and Tvedt (1983). 

6.3 Brittle Fracture Reliability Analysis 

In order to model the singularity at the crack-tip, Bestfield et al (1990) used enriched 
element (Gifford and Hilton, 1978) . Other methods such as J-integral approach (Rice, 1968) and 
hybrid elements (Akin, 1976; Barsoum, 1976; Henshell and Shaw, 1975; Tong et al., 1973) can 
also be used. The enriched element approach has the advantage that mode I and II stress intensity 

factors, k t and k are directly calculated along with the nodal displacement. This simplifies the 

development of the sensitivity equations which are needed in first-order reliability analysis. 
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The discretized global finite element equations are obtained by assembling the enriched 
elements that surround the crack-tip and the regular elements that model the remainder of t e 
continuum. The global system of [neq + 2] equations (i.e., number of displacement equations plus 
mode I and II stress intensity factors) is 


K(b) 6(b) = f(b) 


( 6 - 8 ) 


where the generalized displacement, 6, and external force, f, vectors are 



respectively and the global stiffness matrix, K(b), is given by 

R(b) C(b) ' 

C(b) T E(b) . 


K(b) = 


(6-9) 


( 6 - 10 ) 


In Eqs. (6-8) through (6-10): d and h are the regular displacement and external force vectors, 
respectively; R, E, and C are the [neq x neq] regular stiffness matrix, the [2 x 2] stiffness matrix 
from the enriched terms, and the [neq x 2] coupled stiffness matrix from the regular and enriched 
terms, respectively. The other submatrices in Eq. (6-9) are 


K(b) = 


Kj(b) 

K „< b) i 


I f i (b> 1 

and fl(b) = j f ^ (b) j 


( 6 - 11 ) 


where the two terms fj and f }1 are zero if the enriched element is not on a loaded boundary. 
Equations (6-8) through (6-11) are solved by condensing out the stress intensity factors (i.e., static 
condensation). 


For mixed Mode I and Mode II fracture, several kinds of fracture criteria have been 
summarized by Wu and Li (1989). Among these criteria, the most widely used are: the maximum 
principal stress criterion proposed by Erdogan and Sih (1963) and the minimum strain energy 
density criterion, Sih (1974). In the case of mixed mode fatigue, the fatigue laws .are generally 
based on an equivalent Mode I case to simulate actural mixed mode behavior. In order to be 
consistent with the mixed mode fatigue laws, the maximum principal stress cntenon (Erdogan and 
Sih, 1963) is applied here to determine the equivalent Mode I stress intensity factor. Thus, the 
performance function for the mixed mode fracture can be expressed as 


g(b) - Kc * *eq (6-12) 

Equation (6-12) implies that fracture occurs when the equivalent Mode I stress intensity factor, 
Keq, exceeds the critical value, l^. The direction of crack growth where the hoop stress becomes 
maximum is given by 


T 

Z(K, 0) = 0 


K = 0 


(6-13) 
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where 


0 = 


sin 0 

3 cos 0 - 1 


(6-14) 


In Eq.(6-14), 0 is measured from the current crack line. The relation between the equivalent Mode 
I stress intensity factor (k^ ) and stress intensity factor (tq, K n ) is given by 


T 

k = O K 
eq 


where 


O = cos 




3 sin j j 


(6-15) 


(6-16) 


and 0 is determined by Eq. (6-13). When only Mode I or Mode II fracture is present, Eq. (6-12) 
can be rewritten as 

g(b) =k c -Kj , i = I, II (6-17) 


where Kc is given by 


iCc=Kj c (for Model) and Kc-^pKlc (for Mode II) (6 18) 

In Eq. (6-18), kj c stands for the fracture toughness. As indicated in Section (6.2), the 
determination of the reliability index for calculating the first-order probability of failure in the 
FORM is achieved by solving an optimization problem with one constraint (limit-state condition). 
In order to incorporate other constraints such as equation of equilibrium, crack direction law (in 
fatigue crack growth problem) in the formulation, the method of Lagrange multipliers can be 
applied. The statement of the optimization problem for brittle fracture is described in the 
following. 

The nonlinear programming problem consists of determining the correlated random 
T T T T ... 

variables, b = [b , ..., b ] , and the generalized displacements, 5 = [d , k ], that minimize the 

distance from the origin to the limit-state surface in the independent standard normal space. The 
minimizer is termed as the reliability index (3 (Eq. (6-7)). The minimization is subject to the 
following equality constraint: 

K(b) 5(b) = f(b) (6-19) 

(i.e., equilibrium) and the following inequality constraint: 
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( 6 - 20 ) 


g(b) < 0 

(i.e., the performance function being on the limit-state surface is a constraint in the optimization 

Proble^quations (6-19), and (6-20) are converted to the Kuhn-Tucker problem (Arora 1989) by 
defining a Lagrange functional, L, of independent variables b, 8, p, X, and cx as follows 


L (b, 8, p, X, a) = r T r + p [f - K 5] + X [g + a ] 


( 6 - 21 ) 


where p is a Lagrange multiplier for equilibrium, X > 0 is a Lagrange multiplier for the inequality 
constraint, and a is a slack variable that is introduced to ensure that g < 0. Depending on the sign 
of X, the function to be minimized will increase or decrease with a change in g. In other words, if 

X > 0, then r T r will decrease (i.e., minimize) while g < 0 (Converse 1970). The Kuhn-Tucker 
necessary conditions for the minimization of Eq. (6-21) are obtained by setting the derivatives of 

the Lagrange function with respect to the independent variables b, 8, p, X, and a to zero, i.e., 


dL 3[rTr 1 + p T It - K 8]} + A. = 0 


5b 


3b 


( 6 - 22 ) 


— = -p T K + X^=0 

35 38 


(6-23) 


dL 

3p 


= f - K 8 = 0 


* = g + a 2 , 0 

dX 


(6-24) 

(6-25) 


— = 2 X a = 0 
9a 


(6-26) 


The optimization requires the solutions of Eqs. (6-22)-(6-26) for b, 8, jx, A. > 0, and a. Equation 
(6-24) is simply equilibrium; and Eqs. (6-25) and (6-26) can be simplified to eliminate the slack 

variable, a, such that, X g = 0 and g < 0 which ensures that X > 0. 

Since 8 and b are independent variables in the Lagrange function (see Eq. (6-21)), the 
partial derivative of the second term with respect to b in Eq. (6-22) can be expressed as 




(6-27) 


To simplify the right hand side of Eq. (6-27), the first-order probabilistic finite element equation 
(Eq. (3-4)) is employed for the present static problem, i.e.. 


*r 38 _ 3f 3K g 

K 3b-5b-'3b i> 


(6-28) 
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which yields 


3S[f-KS>K§ (6 ' 29) 

when Eq. (6-28) is substituted into Eq. (6-27). Now multiplying each side of Eq. (6-29) by p. T 
and using Eq. (6-23) in the right hand side we obtain 

<6 - 30) 

which can be expressed as 

H T ;|[f-KS]^|i§K <«') 

since g is only a function of k. 

Substituting Eq. (6-31) into Eq. (6-22), the final optimization problem becomes 

L P + X L S = 0 (6-32) 


when X > 0 where 


(6-33) 


(6-34) 


r) T 

In Eqs. (6-34), ^ [r r] is computed either explicitly or by finite difference depending if the 
random variables are normal or non-normal. In order to perform the sensitivity analysis on the 

stress intensity factors, namely, the probabilistic finite element method described in Sec. 3 can 

be applied to accomplish this task. Since we are only interested in the sensitivity of the stress 
intensity factors, considerable computational effort can be saved by using the adjoint method as 
described in Section 5.2. The iteration algorithm for the brittle fracture reliability is given by 
Besterfield et al. (1990). 

In order to demonstrate the applicability of this approach to the brittle fracture reliability 
analysis, a single edge-cracked beam subjected to a concentrated point load is considered (see Fig. 
7-5). The problem constants are given in Table 7-1. Due to symmetry, ten regular 9-node 
elements and two enriched 9-node elements are depicted in the left half of the beam as shown in 
Fig. 7-5. The applied load is modeled with one random variable with a coefficient of variation of 
0.1 and the crack length is also modeled with one random variable and the coefficient of variation 
of 0.1. The convergence criterion for the optimization is 0.001. The variance of the Mode I stress 
intensity factor with randomness in force, material, crack length and the combination is presented 
in Table 7-2 for the adjoint method. Also presented in Table 7-2 are the summaries of the 
numerical performance and results of the reliability analysis (e.g., starting point, number of 


■ g _ . dg. 3 k 

L ‘ 3b + 9ic 3b 


, P 3 , t . 

L = 3i: [r r] 
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iterations, the failure point, reliability index, and probability of failure). As shown in Table 7-2, 
the random crack length has less effect on the probability of failure due to the smaller variance for 
random crack length. 

6.4 Fatigue Crack Growth Reliability Analysis 

Fatigue crack growth is sensitive to many parameters and these parameters can seldom be 
determined accurately. Uncertainties in the crack geometry, material properties, crack direction, 
crack length, component geometry, and load time history all play a role. Thus, the prediction of 
fatigue failure must be treated as a probabilistic problem. 

The first order second moment reliability method (FORM) can be applied to this problem as 
before by solving a constrained optimization problem. Due to the combined effects of external 
loading, unsymmetrical component geometry and crack geometry, cracks rarely grow in a straight 
line. Thus, the mixed-mode fatigue crack growth law and crack direction law should be employed. 

The most common law for fatigue crack growth is the Paris-Erdogan model (1963), which 
gives the fatigue life, T, by 


T r da (6-35) 

. .n 

D (Ak ) 

J eq 

a^ 

where a\ and af are the initial and final crack lengths, respectively, da is the random crack path; D 
and n are primarily material parameters but can also depend on the loading and environmental 
effects; and Ak (a) is the range of equivalent Mode I stress intensity factors, i.e., 


max min 
Ak = k - k 
eq eq eq 


(6-36) 


where k and k are the minimum and maximum equivalent Mode I stress intensity factors, 
eq eq . . 

respectively, associated with the minimum and maximum cyclic applied stresses, respectively. Ir 
the minimum equivalent Mode I stress intensity factor is assumed to be zero, then 


max 

Ak = k 
eq eq 


K 

eq 


(6-37) 


The direction of the crack can be considered to be a random function, which will depend on the 
material properties and the history of the loading and the crack path. At each step, the statistics of 
the crack-tip, as reflected in this random function, in conjunction with the previous length of the 
crack and its orientation, will be used to obtain the new configuration. Based on the maximum 

hoop stress criterion (Erdogan and Sih, 1963), the crack growth direction Z(k, 0) given by Eq. (6- 
13) is also employed here. Unlike the case of brittle fracture discussed before, the performance 
function for fatigue crack growth is given by 

g = T - T (6-38) 

5 s 
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where T is the service life of the component. In other words, the component fails when the 

s 

fatigue life is less than the desired service life. The performance function could also be expressed 
in terms of a critical crack length. 


The calculation of the reliability index by the first-order probability theory is performed in 
the same way as before by solving a constrained optimization problem. Before stating the 
optimization problem for fatigue crack growth, the crack direction law (Eq. (6-13)) must be 
discretized into "npts" discretization points along the crack path. At each crack path discretization 
point, the crack direction is 

Z. = k, — 0 k = 1, -, npts (6-39) 

k k k 


where and represent k and © evaluated at £ = k = 1, ..., npts. Thus, at each crack path 

discretization point, the new crack direction is recalculated and the crack is then allowed to grow to 
the next discretization point. 

The calculation of the reliability index by the ’first-order probability theory is posed as a 
constrained optimization problem. Unlike the previous brittle fracture reliability problem, both 
equality and inequality constrains have to be satisfied at each crack path discretization point. Also, 
the crack direction law (Eq. (6-39)) has to be included in the Lagrange function. By defining a 

Lagrange function, L, of independent variables b, Pj, ..., P- npts > ..., $ n p ts , ’ ^npts’ 

0,, ..., 0 , X, and a, we have 

1 npts 

L (b, p., 8., <p., 0., X, a) 


T 

= r r + 


npts 

X M f i 
1 = 1 


K. 8. ] + 
i i 


npts 

X* 

k=l 


k k 


+ X T - T + a 


(6-40) 


where p. is a Lagrange multiplier for equilibrium, (p^ is a Lagrange multiplier for the crack 

direction law, a is a slack variable which is introduced to ensure that g < 0 and X > 0 is a Lagrange 
multiplier for the inequality constraint. 

The necessary conditions of Eq. (6-40) (i.e., derivatives with respect to the independent 
variables) then lead to 


3L 

3b 


npts npts 

^tr T r] - X 1*1 1 l f i - K i 5 il + X Wb +X T 'b " # 

i=l k=l 


(6-41) 
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9L 



f. - K. 8. 

= 0 


i = 1. 

..., npts. 

no sum on i 

(6-42) 

^i 


i i i 










npts 






dL 

98. 

1 

= 

T 

- p. K. + 

M 1 

t 

k=l 

*k ®k ) 'S. + 
1 

XT -8. 

1 

= 0 








i = 1, 

..., npts. 

no sum on i 

(6-43) 

dL 


2 

T-T + a 

= 0 





(6-44) 

ax 


s 







dL 


2 X a = 0 






(6-45) 

da 









dL_ 


Z. = 0 




i = 1, 

..., npts 

(6-46) 

9(p. 


i 









npts 

9K. 

npts 

• 




9L 

90. 

1 

= 

k=l 

k 

90. 

i 

8 k - % ^k 
k=l 

( \h. 

i 

+ X T, 0 

1 

0 








i = 1* 

..., npts 

(6-47) 

where: T,. 

b 

in 

Eq. (6-41), 

T, g in Eq. (6-45), and T, 0 in Eq. (6-47) 

are the derivatives of the 


fatigue life T with respect to b, 8., and 0., respectively, assuming b, 5., and 0. are independent 
variables; and (Z R ), b in Eq. (6-41), (Z^g in Eq. (6-43), and (Z k ), 0 in Eq. (6-47) are the 

derivatives of the crack direction law Z k with respect to b, 8., and 0., respectively, assuming b. 


8., and 0. are independent variables. The optimization requires the solutions of Eqs. (6-41)-(6-47) 

for b. H . 8, 8 npts , a.XiO.9, 9„ pts and .... 8^. Equation (6-42) 

is simply equilibrium at each discretization point; Eq. (6-46) is the crack direction law at each 
discretization point; and Eqs. (6-44) and (6-45) can be simplified to eliminate the slack variable, ct, 
such that X g = 0 and g < 0, which ensures that X > 0. 


Since 8. and b are independent variables in the Lagrange function [see Eq. (6-40)], the 

partial derivative of equilibrium with respect to the correlated random variables in Eq. (6-41) can be 
expressed as 

9f. dK. 

^ [f. - K. 8.] = ^ 8. no sum on i, i = 1, .... npts (6-48) 

To obtain an expression for the right hand side of Eq. (6-48), the first-order probabilistic finite 
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element equation (see Section 3) is employed, i.e., 

35. 3f. 3K. 

K — l - = — l - - — - 5 

i 3b 3b 3b °i 


no sum on i, i = 1, npts 


(6-49) 


which yields 

j3_ 

3b 


35. 


[ f , - K i 5 ,1 ■ K i 3b 


no sum on i, i = 1, ..., npts (6-50) 


when substituted into Eq. (6-48). Since 5^, ..., 8 and, 0^, ..., 6^^ are independent 
variables in the Lagrange function, Eqs. (6-43) and (6-47) simplify to: 


- Hi K i + <Pj ( z j)>5 + ^ T ’§ = 0 


no sum on i, i = 1, ..., npts (6-51) 


and 


3K. 


- u. — - 5. + (p. (Z.), n + XT, q =0 no sum on i, i = 1, ..., npts (6-52) 

i 36. 1 1 1 


respectively. Multiplying each side of Eq. (6-50) by |i. and substituting in Eqs. (6-51)-(6-52) 
yields 


r 3K. 

T ~ K'. 1 — -5. - T 


■\ 


>*[ k [f , - K , 8 i ] ■ 


i 38. i ■»; 

3K + T ’ 


- u a'w , 

(Z.), 0 -(Z.), R K'. 1 —5. 
L’6 i 5. i ae. i 


5. 


35. 
^ 3F 


no sum on i, i = 1, ..., npts (6-53) 


which can be expressed as 


*1 & i f i - K i 8 i ] 




3K. 

T >s K — 5 - - T, 0 
’ 8 : 1 36 . 1 e ; 

1 i 1 

3K. 

(Z.), 0 -(Z.), s K". —5. 

^ i A e. i A 5. i ae. 1 


(z.), + t, 

v r’K. ’k. 

i i 


3k^ 
? 3b 


no sum on i, i = 1, ..., npts (6-54) 
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since T and Z. are only functions of k., i = 1, .... npts. After substituting Eqs. (6-52) and (6-54) 
into Eq. (6-41), the final form of the optimization for reliability for fatigue crack growth is given 
by 


. P - . g _ n 

L + A. L — 0 


(6-55) 


where 




npts 

-I 

i=l 


. 9K. 

T.~ K'. 1 — 8. - T 


’S. i ae. i *e 
1 3R— < z s >- + T ’ 


> 


u . 

(Zi), 9 -(Zi), 8 K-. ^5. 
1 1 1 


r 7 k. 


K. 

i 


9k. 

i 

9b 


npts 


1 9K. 

w l < 2 i> e. - < z i>-6. K 'i ^ 8 i 

i i i 




> + T. 


(6-56) 


(6-57) 


to Eq. (6-57): T T , T >e , T.g , (Z.) v (Z,)^. (Z k ), 0 , and (Z.), 5 are determined 


1 1 


explicitly in Reference (Besterfield et al., 1991). The sensitivity of the stress intensity factors, 
3k. 

namely , is also computed by the PFEM. 


In order to demonstrate the performance of the method for reliability analysis against failure 
due to fatigue crack growth, a classical Mode I fatigue problem is presented. Figure 7-6 shows a 
finite rectangular plate with a single edge crack of length a subjected to a distributed load. The 
problem constants and second-moment statistics are given in Table 7-3. Due to symmetry, two 
enriched 9-node elements and twenty-three regular 9-node elements are depicted on the upper half 
of the plate. The reliability index is plotted versus the service life under the various types of 
uncertainties for the reference solution 13 and the solution obtained by PFEM in Figs. 7-7a and 7- 
7b, respectively. The same trends as the reference solution with the slight difference in the value 
of the reliability index can be observed through comparison of Fig. 7-7a with Fig.7-7b. This 
difference is due to the small numerical error in calculating the stress intensity factor by finite 

element methods. As shown in Fig.7-7a, for a service life of 4 X 10 6 cycles, the reliability index 
is less for uncertainty in the initial crack length (100 % coefficient of variation) and stress (25 % 
coefficient of variation ) than for randomness in the final crack length (10 % coefficient of 
variation), fatigue parameter D (30 % coefficient of variation), and fatigue parameter n (2.5 % 
variation). When all five of the parameters are treated as random, the combined effect is much 
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greater than any one individual effect, as expected. 

7. SBEM FOR THE CURVILINEAR FATIGUE CRACK GROWTH RELIABILITY ANALYSIS 
7.1 Introduction 

The development of probabilistic finite element method (PFEM) and its applications to 
linear, nonlinear structural mechanics problems and fracture mechanics problems have been 
discussed in the previous sections. In this section, we present a novel computational tool, called 
Stochastic Boundary Element Method (SBEM), for the reliability analysis of a curvilinear fatigue 
crack growth. 

The SBIM based on the perturbation techniques has been developed by Ettouney et al. 
(1989) and Dasgupta (1992) for quantifying the statistic moments of tractions and displacements of 
a stochastic system. A general methodology, which combines the first order reliability method 
(FORM) with the mixed boundary integral equation method (Lua et al. , 1992c), has been 
formulated most recently by the authors (Lua et al., 1992d). The preformance and efficiency of the 
developed SBEM have been demonstrated by the problem of probabilistic fatigue crack growth. 

The state-of-the-art of boundary element methods along with various of codes are given in 
the Boundary Element Reference Book (Mackerle and Brebbia, 1988). Due to its modelling 
efficiency and solution accuracy, BEMs have been used extensively in the field of computational 
fracture mechanics (Aliabadi and Rooke, 1991; Cruse, 1988). The application of the BEM to a 
curvilinear fatigue crack growth is presented in this section. 

The curvilinear fatigue crack path is mainly attributed to the inherent inhomogeneity of the 
advance materials such as ceramics, composites or polycrystalline alloys. The existence of a micro- 
defect such as a void, a rigid inclusion or a transformation inclusion perturbs the stress field at a 
growing crack tip, resulting in a curvilinear crack path. In order to model the singularity at a moving 
crack tip, an automatic remeshing in conjunction with the quarter-point singular element (Barsoum, 
1976; Henshall and Shaw, 1975 ) has been developed by Saouma (1984) to study the fatigue life of 
attachment lugs. A remeshing scheme based on the Arbitrary Lagrangian Eulerian (ALE) together 
with enriched finite elements has been developed by Besterfield (1991) in the reliability analysis of a 
fatigue crack growth. For problems of multiple fatigue cracks in which elastic interactions of a 
fatigue crack with micro-defects are treated, the remeshing scheme will be prohibitively complicated. 
The formulation based on the Boundary Integral Equations (BEEs) has several advantages in terms 
of solution accuracy and modeling efficiency. 

Due to the degeneration of the usual displacement BIE for coplanar crack surfaces, the 
traction BEE has to be employed on the crack surface. The traction BIE alone is insufficient to solve 
the problem due to the coupling and interaction of the boundary of the component with the growing 
crack. Thus, the displacement BIEs have also to be applied. This set of mixed BIEs provide a 
unique solution for the boundary value problem. The application of the mixed BIEs to the elastic 
interactions of a fatigue crack and a micro-defect can be found in the Reference (Lua et al., 1992c). 

By adding a few elements to permit crack extension along the crack growth direction, 
remeshing can almost be avoided. Similar to the approach used in enriched finite elements (Gifford 
and Hilton, 1987), a special interpolation function which incorporates the stress intensity factors is 
employed to model the near tip Crack Opening Displacements (CODs). The mixed BIEs are 
presented in Sec.7.2 for a multi-connected region with a fatigue crack. An enriched element which 
incorporates the mixed mode stress intensity factors is applied to characterize the singularity at a 
moving crack tip. The response gradient, which is key in FORM, is determined in Sec. 7.3 by 
direct differentiation. Due to the presence of three random processes in the expression of the 
response gradient, namely the mode I and mode II stress intensity factors and the crack direction 
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ancle the first order response- surface model is employed to determine the response sensitivity of 
these random processes. An iteration scheme based on the HL-RF method (Rackwitz and Pass er 
1978) is employed to find the most probable failure point (or design point). Due to the high 
accuracy of the response gradient calculation based on the direct differentiation, fast convergence is 
obtained in the numerical iteration. The accuracy and efficiency of the present approach are 
demonstrated in Sec. 7.4 through a fatigue crack growth problem with randomness in the crack 
geometry, defect geometry, fatigue parameters and external loads. 

7.2 Mixed BIEs for a Multi-Connected Region 


Figure 7-8 shows a finite linear elastic body bounded by outer boundary T 0 and inner 

boundaries T; (i = 1, 2, ..., M), containing a finite crack under remote loading t*. A local 

Cartesian coordinate system (x',y) with origin at the center of the crack is employed with the y - 

axis normal to the crack surface r c . On the displacement boundary T u , the displacement 

components u* are prescribed; and on the remaining outer boundary T t , the traction components t, 

M _ 

are given. The boundary conditions on the inner boundary r, = ^ Tj can be specified based on 

i=l 

the characteristics of a micro-defect (Lua et al., 1992c). 

The usual displacement BIEs (see e.g. Telles, 1983) can be successfully applied on both 
inner boundary H and outer boundary r o . The resulting BIEs on H u T 0 are 


Cij Uj(Q = 


u- k (C;x) t k (x) d f(x) - f 4(C;x) u k (x) dT(x) 


n*r 0 


n-r 0 


Tin j C k (C ;x') Au k (x') dT(x') 
for^Gr,ur o (7-1) 


where the symbol "f stands for the principal value of the integral in Cauchy s sense, and all the 
quantities with the prime ( ' ) are defined in the local coordinate system; t k (x) and u k (x) are the 
components of traction and displacement, respectively, in the global coordinate system. u ik(C; x ) 
and 4(C;x) represent the displacement and traction, respectively, in the k-th direction at the field 
point x corresponding to a unit point force applied at the source point £ in the i-th direction. 
Explicit expressions for these free space Green’s functions are given in Telles (1983). The 
transformation from global to local coordinate system is given by the following transformation 

matrix T (or T j n ): 


T = 


cos0o - sinBo 


(7-2) 


L sinQo cosQo J 

The coefficient matrix Cjj depends on the smoothness of the boundary; = £ 5jj (for smooth 


boundary). The quantity Au m (x') designates the COD on the crack surface T c in the local 
coordinate system defined by 
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Au m (x') = Uf^(x') - Ufn(x') 


(7-3) 


where u^(x') and Ufn(x') are the components of the displacement on the upper and lower crack 
surface, respectively. The coupling term representing the effect of CODs, which differs from the 
usual displacement BEs, is the 3rd term in Eq. (7-1). 

The displacement BIE (7-1) alone is insufficient to solve for all the unknowns, namely, the 

unknown displacement and traction on H u T 0 and CODs on T c . Due to the degeneration of the 
displacement BIEs for coplanar crack surfaces, the higher order BIEs based on the traction 

boundary data on T c is employed. Using the prescribed traction boundary condition on T c , the 
resulting traction BEs become 

ti(0 = nj { | u* jk (C;x) t k (x) dT(x) - v!jk(C;x) Uk(x) dHx) - 

• n ♦ n o 

. Y Vnk(C **) Au k (x) dr(x ^ 5 for C e r c ( 7 ' 4 > 

J r c 

where the symbol ( f) stands for the finite-part of a divergent integral, t; are the prescribed traction 
components on the upper crack surface T + c , and nj=(sin0 o , -cos0 o ) are the components of the 
normal vector to V + c in the global coordinate system. The free space Green's function u tjk(C;*) and 
Vijk^C>x) are given by Lua et al (1992c). 

In order to characterize the crack tip singularity, an enriched element which incorporates 
the stress intensity factors (SEs) is used at the crack tip: 


' n-r„ 


* Tj n T JS 




(7-5) 


where s is the distance behind the crack tip, a is the semi-crack-length, and Ki and K n are mode I 
and mode II SEs. In the numerical implementation of the mixed BIEs (Eqs. (7-1) and ( 7 - 4 )), all 

the boundaries, namely Tj , T 0 and T c have to be discretized first. By dividing the boundary 
Ti + T 0 and the crack surface T c into NE and NC elements, respectively, the discretized version of 
the mixed BEs can be expressed as 


NE r t NC 

Cij Uj(0 - I [GSKt;*m) I? - H: k (C;x m ) u|»J - 1 QSE(C ;*»> au 


m 

k 


m=l 


ti(C) = nj(C) | X, [Dijk(C; x m) Ik ■ Sy k (C; x m) u k*] ' 2, ^ijk^C > x m ) Au k 

m=l 


m=l 

for C, = xi, X 2 , ..., xne 
NC 


(7-6) 


m=l 


for £ = X! , x 2 , ..., x N c 


(7-7) 
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where the coefficient matrices G[" k , Fft , Q* , D™ k , , and are given by Lua et al (1992c). 

7.3 FORM for the Curvilinear Fatigue Crack Growth 

As described in Sec. 6.2, for the general correlated and non-normal random variables b, 
three steps are required in the first-order reliability analysis (FORM). They are: 1) transformation 
of b into the uncorrelated standard normal vector r by Rosenblatt transformation (Rosenblatt, 
1952), 2) approximation of the failure surface in the r-space by a flat hyperplane at the most likely 
failure point (design point), 3) determination of the reliability index b by computing of the 
minimum distance from the origin to the limit state surface. As discussed in Sections 6.3-6.4, the 
design point has to be determined by an iterative optimization algorithm, which involves repeated 
computation of the limit state function (Eq. (6-1)) and its gradient, in order to ensure rapid 
convergence, an accurate determination of the response gradient is required. 

Unlike the previous sections (6.3 and 6.4), where the response gradients or sensitivities 
have been determined by PFEM, the response gradients are calculated using SBEM. The direct 
differentiation coupled with the response-surface method is employed to perform the sensitivity 
analysis. 

Assuming that the crack geometry (a; , a f , x 0 , yo , 9o)> fatigue parameters (D, n), external 
load (x), and defect geometry (x c , y c , r c , pj) (see Fig. 8) are modeled by a q-dimensional random 
vector b, the performance function for a fatigue problem is given by Eq. (6-38). Since the service 
life T s is a deterministic variable, the gradient of the limit state function is given by the response 

sensitivity, i.e., 


ff=& T < b -'W>» < 7 - 8 > 

where T is given by Eq. (6-35). In order to facilitate the response gradient calculation, the line 

mapping is applied to map a curvilinear crack path to a local coordinate system, ^ ( ^ e [-1 , +1] ). 
The mapping function is defined by 


a = y[( a f- a i)4 + ( a f + a i)] (7-9) 

Assuming that the crack geometry (a; , af , xo , yo , ©o) > fatigue parameters (D, n), external load (x 
), and defect geometry (x c , y c , r c , pi) (see Fig. 8) are modeled by a q-dimensional random vector 
b (or r in the transformed space), and using Eqs. (6-37) and (7-9), Eq. (6-35) can be rewritten as 


T(b, K eq (b)) = £’f(b, K^b, 5)) d^ 
where the function f in Eq. (7-10) is given by 


f(b, K eq (b,^) = 


J(b) 


D[Keq( *»,$)]" 

The Jacobian J in Eq. (7-11) is defined by 


(7-10) 


(7-11) 


J(b) = £(a f -ai) 


(7-12) 
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Note that the function k (b, £) can be described only in an algorithmic form through the mixed 
BIEM (Eqs (7-1) and (7-4)). This is the place where the mixed BIEM interfaces with the FORM 
to form the stochastic BIEM. Using Eqs. (7-10), (7-1 1) and (6-15), the total denvauve of the 

8T 

response ( — ) is given by 

9b 


8T _ r 

<*> u 


f - b + f > K «<, * K «u, 0b 


atc^b, ^) 


+ K 


8K u (b, 4) 80(b, 5) 

+ K, 


«U n 8b 


«t.e 8b 


(7-13) 




where K„„ , and K„ . are derived from Eq. (6-15) and given by 

3 9 

K «U, = COS 2 


«t.K u 


~ 2 e - e 

= -3 cos — sin — 


3 e 

Ve = 2 C ° S 2 


0.0 L . 2 0 2 ®] 

■ cos- sin- Kj + 12 sin y-cos 


(7-14) 

(7-15) 

(7-16) 


Both f b and f K in Eq. (7-13) can be determined explicitly from Eq. (7-11). The results are 


f,b = 


1 9J 


8D J ^(K^) 8n 


f K =- 


D(kJ" 3b ' D 2 !^" 3b ' D (k,J" 3b 


(7-17) 


(7-18) 


D 


Due to both the complicated explicit expressions and implicit functions involved in Eq. (7- 
13), numerical integration is required to calculate the response sensitivity (Eq.(7-13)). By dividing 
the integration interval [-1, +1] into Npts-1 line elements, which correspond to Npts-1 crack 
growth steps, and applying the trapezoidal rule, Eq. (7-13) can be approximated by 


3T JS? 


f .b + f .K„ ( Vk, 


8Kj(b, 4) 

8b 


+ K 


8K u (b, ^) 

8b 


+ K 


30(b£) 


8b 


(7-19) 


W„ 


where = -1 + 2(n-l)/(Npts-l) , and W n are the integration weights given by 

W n = — -i (for n=l or n=Npts) , W n = — 2 — (otherwise) 

Npts-1 Npts-1 


(7-20) 


Since Kj, K n and 0 are implicit functions of both b and % in Eq. (7-19), the direct calculation of 


31 



3kj aK„ 30 . , ., . . ,. 

these response sensitivities, namely, and ^ is not feasible m this case. 

As shown in Eq. (7-19), the key step in the implementation of stochastic BIEM for a 
curvilinear fatigue crack growth reliability is to determine three response sensitivities ( 

— 1 __ii and ) at each integration point (n= 1, 2, .... Npts). In addition to the implicit 
3b ’ 3b 3b 

dependence of functions Kj(b, ^), iCjj(b, and 9(b, on both b and these three functions 
represent three random processes due to the presence of the random vector b. Since 
Kj(^), Kjj^) and 0(£) at a given realization of b can be easily generated by using the mixed BIEM 
(see Sec. 7.2), the response-surface approach (Faravelli (1989, 1986)) is employed to determine 

the response sensitivity ^ and at each discretization point (n =1 , 2 Npts). The 

first order response model in b is employed in conjunction with factorial experiments with each 
factor at two levels (Myers, 1971). As K,(b, ^), K n (b, ^) and 0(b, ^) are independent of the fatigue 

parameters, D and n, the dimension of b in Kj, Kq and 0 is q-2. 

With the help of the deterministic solver based on the mixed BIEM, 2 q ' 2 computer 
simulations are performed in accordance with predefined factorial simulations. By performing the 

least squares fitting process at each crack path discretization point (n =1, 2, ..., Npts), the 

3Kj(b, ^) (a^b,^) 30(b 

history of the response sensitivity of — ^ f 

be determined. 


and 


can 


at the L-th iteration, b L 
the results 


3b 3b 

Substituting 

(3 K i(b’ jjj, ^ and ; n = 1, 2, ..., Npts) into Eq. (7-19), the sensitivity of the 

3b ’ 3b 9b 

fatigue life T with respect to the primary random vector b at b L can be determined. For a given 
service life T s , an iterative algorithm to obtain the location of the design point b , the response 

[ 9t ] 

and the reliability index (or the probability of failure) can be 


sensitivity at the design point 


3b 


found in the Reference (Lua et al., 1992d). 


7.3 Numerical Results 

In order to show the accuracy and efficiency of the stochastic BIEM in a curvilinear fatigue 
crack reliability setting, a single edge-cracked plate with a random transformation inclusion is 
considered (see Fig. 7-9). The plate geometry (W, L), initial crack location (x 0 , yo), initial crack 

angle (0o), final crack size (af), and material constants (aluminum 7075-T651) are deterministic 
parameters given by 

L = W = 2.0 in ,xo = -1.0 in ,yo = 0.0 in , 0o = 0.0 (7-21) 

af = 0.5 in , p. = 3.866 x 10 6 psi , v = 0.33 (7-22) 

where is the shear modulus and v is the Poisson's ratio. The crack geometry (ai ), external load 
(t), fatigue parameters (D, n), the defect geometry (x c , y c , r c ), and the internal pressure (Pi) 
resulting from the residual strain in the inclusion are assumed to be independent random variables 
with specified probability density functions. 

The statistical parameters of random input variables (mean, standard deviation and 
coefficient of variance (COV)) along with corresponding distribution functions are listed in Table 
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1. As shown in Table 4, the initial crack size aj has the largest dispersion (COV =60%). For the 
initial crack length a;, a uniform distribution with a tail is employed here (see Fig. 7-10). The 

detection threshold, which is equal to 7.5 x 10' 3 shown in Fig. 7-10, represents the lower limit of 
an inspection device to detect the presence of a small crack. Below the detection threshold the 
probability density is assumed uniform; above the threshold the probability density decays linearly 
to zero, representing false negatives of the inspection technique. For the purpose of verifying the 
accuracy of the stochastic BIEM, the Monte Carlo Simulation (MCS) for the sample size 
N s = 2000 is used. 

The Cumulative Distribution Function (CDF) of the fatigue life T obtained by the stochastic 
BIEM for values of service life T s are presented in Fig. 7-11. The agreement of MCS and SBIEM 
results shown in Fig.7-1 1 demonstrates the accuracy and efficiency of the stochastic BIEM. As a 
rule of thumb (Bjerager, 1989), the sample size necessary for MCS to get a probability estimate 

with good confidence is around 100/pf. For small probabilities of failure pK= 10' 3 - 10' 6 ), which 

are the major interest in reliability engineering, one needs 10 5 - 10 8 Monte Carlo simulations to 
achieve good confidence. The number of iterations in the stochastic BIEM required to find the 

design point b* is only of order 15 to 20 for (3 = 3-5 (or pf = 0.001 - 0.3 x 10" 6 ). Therefore the 
stochastic BIEM based on FORM has an overwhelming advantage over the MCS for small 
probabilities of failure in terms of solution accuracy and efficiency. 

The reliability index (|3) versus the service life (T s ) is shown in Fig. 7-12 along with the 
results of no micro-defect. As shown in Fig. 7-12, the presence of a random transformation 
inclusion has a detrimental effect on the fatigue life. The comparison of response sensitivities at 
Most-Probable-Points (MPPs or design points) versus the probability of failure for both cases is 
plotted in Fig. 7-13. As shown in Fig. 7-13, the presence of a random transformation inclusion 
changes the response sensitivity of a; significantly. The comparison of the loci of the Most- 
Probable- Point (MPP) of crack geometry (a;) is shown in Fig. 7-14. Due to the presence of the 
random transformation inclusion, the locus of MPP of a; changes considerably as shown in Fig. 7- 
14. When the value of a; increases, the probability of failure Pf becomes large (see Fig. 14). This is 
the main reason why the routine crack inspection is so important to avoid the large probability of 
failure. 

8. Conclusions 

The Probabilistic Finite Element Method (PFEM) is presented, which is based on the 
second-order perturbation. Due to the discrete nature of the finite element formulation, the random 
field has to be also discretized. Existing approaches for representation of random fields are outlined. 
To the efficient characterization of the random field, the transformation of the original random 
variables into a set of uncorrelated random variables is introduced through an eigenvalue 
orthogonalization procedure. Both single-field variational principle and three-field Hu-Washizu 
variational principle are employed to develop the PFEM for linear, and nonlinear problems, 
respectively. The computational aspects in the numerical implementation of the PFEM are also 
presented. 

The accuracy and efficiency of PFEM in quantifying the statistic moments of a stochastic 
system are demonstrated through two examples: 1) a stochastic spring-mass system under sinusoidal 
excitation; 2) a cantilever beam subjected to large deflection. The results are in good agreement with 
Monte Carlo simulations (MCSs). The computational efficiency of PFEM far exceeds MCS. Since 
the PFEM developed essentially involve solution of a set of deterministic problem, it is easily 
integrable into any FEM based code. 


33 



The PFEM coupled with the first-order reliability method is alsopresented for the reliability 
analysis of both fracture and fatigue. The methodology consists of calculating the rel^bilny index 
via an optimization procedure, which is used to calculate the probability of fai u ... 
provides^ powerful tool for the sensitivity analysis, which is required m an iterative optimization 
algorithm. Performance of the methodology presented is demonstrated on a single edge-cracked 
beam with a concentrated load and a classical mode I fatigue crack growth problem. 

In addition to the PFEM, the stochastic boundary element method (SBEM), which combines 
the mixed boundary integral equation with the first-order reliability method, is also presented for t e 
c^Xear ^?augu« ra ck reliability analysis. Due .0 the high degree of complexity and nonlmeanty 
of the response, direct differentiation coupled with the response-surface method is employed to 
determine the response gradient. The reliability index and the corresponding probability of failure 
are calculated for a fatigue crack growth problem with randomness m the crack . geometry, < defect 
geometry, fatigue parameters and external loads. The response sensitivity of each primary random 
variable at the design point is also determined to show its role in the fatigue failure. The resu s 
show that the initial crack length is a critical design parameter. Since crack lengths ^ J* 
threshold of an inspection limit are likely to exhibit a large amount of scatter, this makes it imperati e 
that the life expectancy of a strycture be treated from a stochastic viewpoint. 

Probabilistic analysis is becoming increasingly important for the safety and reliability of an 
aging structure and for tailoring new advanced materials. Due to the complexity in ch^actnzing 
material behavior, structural response, and failure mechanism, probabilistic mechanics problems 
computationally intensive and strain the resources of currently available computers. Since many 
sources of parallelism are inherent in probabilistic mechanics problems, it is evident that the 
development of a parallel computing enviroment for probabilistic response analysis in the current 
trend in stochastic computational mechanics. 
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Table 7-1. Problem Constants: Single Edge-Cracked Beam with an Applied Load 


Parameter 

Mean 

Standard Deviation 

Percent 

Length (L) 

10.0 in 

0.0 

0.0 

Width (W) 

5.0 in 

0.0 

0.0 

Thickness (t) 

1.0 in 

0.0 

0.0 

Young's Modulus (E) 

30,000.0 ksi 

3,000.0 

10.0 

Poisson's Ratio (v) 

0.30.0 

0.0 


Applied Load (P) 

10.0 kip 

1.0 kip 

10.0 

Crack Length (a) 

0.1 in 

0.01 in 

10.0 

Stress Intensity Factor (Kj) 

33.452 ksi VS 

0.0 ksi VS 

0.0 

Fracture Toughness (<j c ) 

43.0 ksi VS 

0.0 ksi VS 

0.0 
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Table 7-2. 


Numerical Performance in Brittle Fracture Reliability Analysis 


Randomness 

in 

Variance 

2 3 

(kip /in ) 

Starting 

Point 

Number 

Iterations 

Failure 

Point 

Reliability 

Index 

Probability 
of Failure 

Force 

12.3248 

9.94% 

P=12.5 kip 

5 

P=12.2 kip 

2.173 

1.49% 

Young's 

Modulus 

0.003087 
0.18 % 

mm 

m 

0 

■ 

oo 

0.0% 

Crack Length 

1.8273 
3.83 % 

a=1.4 in 

7. 

a=l-29 in 

2.911 

0.1801 % 

Combined 

i 

12.5107 
10.01 % 

P=12 kip 
E=30e3 ksi 
a=l.l in 

9 

P=12.1 kip 
E=30e3 ksi 
a=1.02 in 

2.079 

1.88 % 
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Table 7-3. 


Problem Constants: Single Edge-Cracked Plate with a Distributed Load 


Parameter 

Mean 

Standard Deviation 

Percent 

Length (L) 

10.0 in 

0.0 

0.0 

Width (W) 

4.0 in 

0.0 

0.0 

Thickness (t) 

1.0 in 

0.0 

0.0 

Young's Modulus (E) 

30,000.0 ksi 

0.0 

0.0 

Poisson's Ratio (v) 

0.3 

0.0 

0.0 

Applied Stress (x) 

12.0 ksi 

3.0 ksi 

25.0 

Initial Crack Length (ap 

0.01 in 

0.01 in 

100.0 

Final crack Length (a^) 

0.1 in 

0.01 in 

10.0 

Fatigue Parameter (D) 

1.0 x 10' 10 

3.0 x 10' U 

30.0 

Fatigue Parameter (n) 

3.25 

0.08 

2.5 
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Table 7-4. The Statistical Parameters 
the Example Problem 


and Distributions of Input Random Variables of 


Random Parameters 

Mean 

Standard deviation 

. cov 

a; (Uniform with tail) 

0.5833 X 10' 2 in 

0.3584 X 10 2 in 

61.4% 

D (Log-normal) 

0.3770 x 10' 9 

0.1885 x 10- 10 

5.0% 

n (Log-normal) 

3.60 

0.18 

5.0% 

x (Normal) 

1 1 .0 ksi 

1.1 ksi 

10.0% 

Xc (Uniform) 

-0.25 in 

0.14433 

57.7% 

Yc (Uniform) 

-0.4 in 

0.05774 

14.4% 

r c (Uniform) 

0.1375 in 

0.03608 

26.2% 

Pi (Normal) 

35.0 ksi 

3.5 

10.0% 
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Figure 7-1 
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Figure 7-3 
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Figure 7-5 
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Figure 7-7a 
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Figure 7-7b 
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Figure 7-12 
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SENSITIVITY OF ai (*E+08) 



Figure 7-13 
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